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Abstract 

We calculate the Hubbard bands for the half-filled Hubbard model on a Bethe lattice with 
infinite coordination number up to and including third order in the inverse Hubbard interaction. 
We employ the Kato-Takahashi perturbation theory to solve the self-consistency equation of the 
Dynamical Mean-Field Theory analytically for the single-impurity Anderson model in multi-chain 
geometry. The weight of the secondary Hubbard sub-bands is of fourth order so that the two- 
chain geometry is sufficient for our study. Even close to the Mott-Hubbard transition, our results 
for the Mott-Hubbard gap agree very well with those from numerical Dynamical Density-Matrix 
Renormalization Group (DDMRG) calculations. The density of states of the lower Hubbard band 
also agrees very well with DDMRG data, apart from a resonance contribution at the upper band 
edge which cannot be reproduced in low-order perturbation theory. 

PACS numbers: 71.10.Fd,71.27.+a,71.30.+h 
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I. INTRODUCTION 



The Dynamical Mean-Field Theory (DMFT) maps lattice models for electrons with a 
Hubbard-type interaction onto effective single- impurity models; for a review, see Ref. 
The parameters of the impurity model must be determined in such a way that the self- 
energy and the Green function of the impurity model agree with the local self-energy and 
the local Green function of the lattice model. The solution of this self-consistency problem 
equally solves the original lattice problem in the limit of infinite dimensions.- For example, 
the Hubbard model on the Bethe lattice with infinite coordination number can be mapped 
onto the single-impurity Anderson model. Then, the self-consistency condition requires that 
its hybridization function and its Green function agree for all frequencies. 

Unfortunately, we are far from an analytical solution of the single-impurity Anderson 
model for a general hybridization function, and a variety of methods have been employed 
to solve the DMFT equations for the single-band Hubbard model. Examples for numerical 
treatments are the Numerical Renormalization Group method, 3 Exact Diagonalization,-"- 
the Random Dispersion Approximation,-^ the Dynamical Density-Matrix Renormalization 
Group (DDMRG) method,-^ and, at finite temperatures, Quantum Monte-Carlo.— - — Ap- 
proximate analytical methods at zero temperature include the Iterated Perturbation The- 
ory^ 3 -, the Local Moment Approach, 14 and the self-energy functional approach.— 

All methods have their merits and limitations and it is desirable to compare their results 
with those from perturbation theory. For the half-filled Hubbard model on a Bethe lattice 
with infinite coordination number, the self-energy^ and the ground-state energy^ are known 
up to and including fourth order in the Hubbard interaction U and up to second order in 
1/C/.— However, these calculations are based on the Hubbard model in infinite dimensions, 
not on the DMFT description. 

In this work, we solve the DMFT equations for the Hubbard model on a Bethe lattice 
with infinite coordination number, Z — > oo, at half band-filling for strong coupling where the 
model describes a paramagnetic Mott-Hubbard insulator. Up to and including third order in 
1/U, we determine the hybridization function of the single- impurity Anderson model which 
corresponds to the Hubbard model on the Z — > oo Bethe lattice. Essential to our approach 
are: (i) the mapping of the single-impurity Anderson model from the 'star geometry' onto 
the 'multi-chain geometry' where each chain represents one of the upper and lower Hubbard 
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sub-bands; (ii) the Kato-Takahashi perturbation theory^^ for degenerate ground states; 
(iii) the Lanczos representation of the hybridization function and the Green function which 
permits an order-by-order solution of the self-consistency equation for the moments of the 
density of states; (iv) the locality of the Hubbard interaction and of the Lanczos operators 
in finite order perturbation theory. 

Our work is organized as follows. In Sect. [Til we introduce the Hubbard model, the 
single- impurity Anderson model, the DMFT equations which link the two models, and the 
two-chain mapping which we use for our perturbative calculations to third order in 1/U. 
In Sect. II I II we adapt the Kato-Takahashi perturbation theory to our problem and use the 
Lanczos algorithm to express the density of states for the (primary) lower Hubbard band 
and for the hybridization function in terms of their moments. Then, the self-consistency 
equation reduces to the condition that the respective moments agree up to trivial signs. In 
Sect. [IV] we investigate the lowest non-trivial order and show how the iterative solution of 
the DMFT equation works in practice. Next, we summarize the results to third order; all 
technical details can be found in Ref . [20] . The remaining problem is the calculation of the 
density of states at the boundary of a semi-infinite chain for a single particle which can 
move between nearest neighbors and experiences a local potential at and near the boundary. 
Its solution and a favorable comparison with previous numerical work-iii is the subject of 
Sect. [V] Conclusions, Sect. IVIj and two appendices, on the secondary Hubbard bands and 
on the Green functions for a particle on a semi-infinite chain, close our presentation. 

II. MOTT-HUBBARD INSULATOR IN DYNAMICAL MEAN-FIELD THEORY 

We start our presentation with the definition of the Hubbard model and the single- 
impurity Anderson model. For a specific choice of the hybridization function in the single- 
impurity Anderson model, its single-particle Green function is identical to the local single- 
particle Green function of the Hubbard model in infinite dimensions. The Dynamical Mean- 
Field Theory prescribes a way to determine the hybridization function self-consistently. 

In general, the single-particle Green function for the single-impurity Anderson model 
cannot be calculated analytically. For the Mott-Hubbard insulator, we use a mapping of 
the model onto a multi-chain geometry where the chains represent the energy levels in the 
energetically separated upper and lower Hubbard sub-bands. 
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A. Hamilton operators and Green functions 



1. Hubbard model 

We consider the repulsive single-band Hubbard model (U > 0) 

H = l if c X^h<y + U ^2 *Wf"a ~ A* ^2 ("*,t + **u) —-T + UD - fiN . (1) 

Here, T denotes the operator for the electron transfer between the lattice sites % and j, the 
fermion operator cf a (cj )Cr ) creates (annihilates) an electron with spin a (=t, i) on lattice 
site i, the operator n i:Cr = c^ CT Cj jCT counts the number of a-electrons on site z, and the operator 
Z) = hi^hi^ counts the number of doubly occupied sites. For the description of the Mott- 
Hubbard insulator, we consider a half-filled system where there is on average one electron per 
lattice site, n = N/L = 1. Moreover, we treat the paramagnetic situation, = = 1/2, 
without any symmetry breaking. The thermodynamic limit, N, L — > oo is implicit in our 
calculations below. 

As a major simplification, we assume that the electrons move between nearest neighbors 
on a Bethe lattice with coordination number Z, 

{—t/yZ if i, j are nearest neighbors , 
(2) 
else . 

Later, we shall let go Z — > oo and choose t — 1 as our unit of energy. The Bethe lattice 
with coordination number Z is an infinite Z-Cayley tree. A Z-Cayley tree is constructed 
from a first site by connecting it to Z new sites which constitute the first shell. One creates 
further shells by adding Z—l new sites to every site in shell s. The Cayley tree has no loops 
and all closed paths are self-retracing.— The Bethe lattice contains s — > oo shells. Since the 
Bethe lattice is bipartite, the chemical potential \i = U/2 guarantees half band-filling at all 
temperatures. 

We are interested in the local Green function of the Hubbard model in its exact ground 
state |^o)- We use the abbreviation 

- W 

for ground-state expectation values and define the local causal Green function in the time 
domain 

G a (i;t) = -i{f B Ci^t)cL (0)}, (4) 



where the Heisenberg operators 

M*) = j*%*e- i&t (5) 
are time-ordered with the help of time-ordering operator T t 

{Q CT (t)ct ,(*') for t > H , 
' { ' h ° 1 ' (6) 
-c+ a ,(t')c ha (t) for t<t>. 

The time-frequency Fourier transformation of the local Green function is defined as 

/OO 
dte^G.fat) (7) 
-oo 

= Inn { (c i)(7 (« - (ff - ^o(iV)) + «?) O 

+ (c+ a (a; + (F- J Eo(iV))-ir / ) _1 Q, CT )} . (8) 



The limit rj — > + is implicitly understood henceforth. In (JS}, E Q (N) denotes the energy 
of the iV-particle ground state |^ ) °f the Hubbard model. The (local) density of states is 
obtained from the imaginary part of the Green function (sgn(x) is the sign function), 

D a (u) = --sgn{u)lm[G lta {uj)] . (9) 

7T 

The density of states is positive semi-definite and its integral over all frequencies is unity.— 
The Green function for non-interacting electrons on a Bethe lattice (U = 0) can be 
calculated in various ways.— In the limit Z — > oo, it approaches the Hubbard semi-ellipse, 



pM = ^wt-(w) 2 for M£H//2 (10) 

with W = 4 as the bare bandwidth. In the presence of interactions (U > 0), the local Green 
function can be expressed with the help of the (proper) self-energy Ti a (co), 



oo 



Note that, in the limit Z — > oo, the self-energy depends only on the frequency- In principle, 
the self-energy can be calculated in diagrammatic perturbation theory.— 
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2. Single-impurity Anderson model 



In order to set up the Dynamical Mean-Field Theory for the half-filled, paramagnetic Hub- 
bard model, we consider the discrete, symmetric single-impurity Anderson model (SIAM) 
in 'star geometry', 

L-2 jj L-2 

#SIAM = ^2 £™ a m,<T«m, CT ~ ^ ^ + 2 ^(<*i<A + ^«m,cr) + Un d ^fl d)l . (12) 
m=0 o" it m=0 cr 

Here, a+ CT (a m ,<r) creates (annihilates) a bath electron with spin a and bath energy £ m , 
dp (da-) creates (annihilates) an electron with spin a on the impurity level with energy 
Ed = —U/2, and fi^ a = d^d a counts the number of cx-electrons on the impurity. The 
Hubbard interaction on the impurity is the same as in the Hubbard model. The parameters 
V m > describe the hybridization between the bath levels and the impurity site. The half- 
filled case corresponds to N = L electrons. The limits N, L — > oo are implicitly understood 
henceforth. 

The single- impurity Anderson model is fully characterized by the hybridization function, 24 

A(w) = — % -— . (13) 

The (causal) Green function of the impurity electrons is defined by 

Gf AM (d;t) = -i(f s l(t)dt(0)) , (14) 

where the expectation value is to be taken in the exact ground state of the single-impurity 
Anderson model. After Fourier transformation, the Green function can be expressed with 
the help of the hybridization function and the (proper) self-energy in the formal 

«r» - jy e ^ m = —^L-— . (15) 

As in the case of the Hubbard model, the self-energy of the single-impurity model can be 
calculated in diagrammatic perturbation theory. 

3. DMFT equations on the Bethe lattice 

The skeleton diagrams for the single-impurity Anderson model and the Hubbard model 



with infinite coordination number are identical; for a review, see Ref 



Therefore, if their 



local Green functions agree, 

GU") = Gf AM (") , (16) 

their self-energies agree as well, 

£» = £f AM (u;). (17) 

The exact solution of the Hubbard model for infinite coordination number reduces to the 
calculation of the Green function of the single- impurity Anderson model for a general hy- 
bridization function A (a;). The DMFT self-consistency equations f|T6|) and (fTT|) single out 
the hybridization function which describes the Hubbard model in the limit of infinite coor- 
dination number. 

For the Hubbard model on the Bethe lattice, the semi-elliptic bare density of states f fTUj) 
results in the following form of the local Green function (llip 

G i;a (z) = ~ (z - V#=t) , z = oj-X a (u) (18) 

so that 

^(u) = u - G L<T {u) - —\- (19) 
holds. Together with f fT5|) . the DMFT equations f|T6|) and f|T7|) reduce to the single condition 

A(u) = Gf AM (u) (20) 

on the hybridization function. The remaining task is to calculate the Green function 
G^ 1 (tu) for the single- impurity Anderson model for a general hybridization function A(co). 
The equation (|2"U|) singles out the hybridization function which describes the Hubbard model 
on the Bethe lattice with infinite coordination number. From now on we shall exclusively 
investigate the single- impurity Anderson model. Therefore, we drop the superscript 'SIAM' 
on all quantities. 

B. Two-chain mapping for the Mott Hubbard insulator 

1. Hubbard bands and charge gap 



We are interested in the description of the Mott-Hubbard insulator where the charge 
gap separates the upper and lower Hubbard bands. Due to particle-hole and spin symmetry 
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(see, for example, Refs. |6lJ25|). it is sufficient to calculate the Green function of the lower 
Hubbard band for a fixed spin, say, o =|, 



For positive frequencies we have Guhb(w > 0) = — Glhb(— oj). Moreover, for the density of 
states we have Z)uhb(^ > 0) = -Dlhb(~ w), i.e., the density of states is symmetric around 
u = 0. 

The upper edge of the lower Hubbard band is the chemical potential \f < for adding 
the Lth electron. The minimal energy for adding another electron to the system (N = L + l), 
the chemical potential is given by fi + = — fi~ .— Therefore, the charge gap obeys 



i.e., it can be obtained from the upper band edge of the lower Hubbard band. 
2. Two-chain single- impurity Anderson model 

The self-consistency equation (|20|) demands that the imaginary part of the hybridization 
function is identical to the density of states. For the Hubbard model at strong coupling, 
U ^> W, we know that the density of states is centered in the two Hubbard bands, \u±U/2\ < 
0(W/2).— For discrete bath levels, the imaginary part of the hybridization function A (a;) 
in eq. ( EEHj) consists of peaks at the bath energies £ m with weights V^. Therefore, the bath 
energies can be grouped into those of the lower Hubbard band, £ m = —0(U/2), and those 
of the upper Hubbard band, £ m = 0(U/2). Consequently, we map the single-impurity 
Anderson model in star geometry, eq. (Tl2|) . onto a two-chain geometry where the impurity 
site hybridizes with two sites which represent the lower and upper Hubbard bands. 28 Note 
that, in numerical treatments of the single-impurity Anderson model, the star geometry is 
usually mapped onto a single chainA Apparently, the two-chain mapping is more adequate 
for the Mott-Hubbard insulator; a similar idea was proposed earlier in Refs. flfl. The 
concept is readily generalized to a multi-chain mapping where each region with a finite 
density of states is represented by its own chain; see below. 

The two-chain mapping can be carried out technically along the lines of the single-chain 




(21) 



A c = jj + - /j, 



2|/i 



(22) 
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mapping (Lanczos tri-diagonalization 3 ). The transformed Hamiltonian reads 

H = H + V = H + Vo + V 1 + V 2 , (23) 

tj (L-3)/2 

^ = ~2 E £(*i!A* - + u <Ml ~ V2)(n d4 - 1/2) , (24) 

Z=0 a 

^° = [(^"0,. + + h.c] , 

(L-3)/2 

Vi = E*i[(«^W + A+A+i,a) + h.c.] , (25) 

(L-3)/2 

^ = X) S e j(«u«i,»- ATA*) ■ 

The d-operators describe the electrons in the lower chain (lower Hubbard band) and the 
/9-operators those in the upper chain (upper Hubbard band). Due to particle-hole symmetry, 
the electron-transfer amplitudes in the lower and upper chain are equal, t~[ = and the 
on-site energies in the lower and upper chains are opposite in sign, e[" = — — e\ — U/2. 
The mapping is shown in Fig. [TJ Later, we shall investigate the model in the strong-coupling 
limit. Therefore, we separated the Hamiltonian into the starting Hamiltonian Ho, eq. fl24|) . 
and the perturbation V, eq. (125|) . Note that H describes the atomic limit, T = 0, where 
there is no transfer between sites in the Hubbard model. 




FIG. 1: Mapping of the discretized SIAM onto two semi-infinite chains, coupled via the impurity. 
The states which have the energy £ m = (U/2) and £ m = —(U/2) in the atomic limit, respectively, 
are mapped onto the upper/lower chain. 



Our task is the calculation of the Green function on the impurity site, eq. ( 121]) . for 
general on-site energies ei and electron-transfer parameters ti. In the two-chain geometry, 
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these parameters assume the role of the energies £ m and the hybridizations V m in the star 
geometry. Our approach relies on the order-by-order expansion of all quantities in 1/U, 



£\ = 

n=0 

oo 

tl = ' 

n=0 



oo , 1 v n 

E'< ' (tj) ■ (*) 

n— n V / 



for I > whereby we implement the self-consistency equation rtzQ 



III. KATO TAKAHASHI PERTURBATION THEORY 

In order to calculate the zero-temperature Green function for the single- impurity Ander- 
son model in strong coupling, we adapt the Kato-Takahashi perturbation theory.— 1 ^ The 
Kato-Takahashi perturbation theory is particularly suitable for Hamiltonians H = H + V 
where the ground state of the unperturbed Hamiltonian Hq is degenerate. 



A. General formalism 



1. Transformation of the eigenstates 

The basic assumption of this perturbation theory is the existence of an isometry r^jy 
from the iV-particle eigenspace £^ of H to the corresponding eigenspace £^ of H, 

Ti,N '■ £i°N ~~ ^ £i,N • (27) 

Explicitly, the Kato-Takahashi projection operator-^"— is given by 

f _ p p(0Vp(0)p p(0)\-l/2 / 98 n 

f it N is unitary provided we may identify the isomorphic subspaces Sf^ and £^n- 

The operators Pf^ project onto the iV-particle eigenstates of Hq, i-e., states with energies 
= E^jf+iU. The operators P^n project onto the iV-particle eigenstates of H. They can 
be expressed as a perturbation series in terms of the projectors Pj^ and the perturbation V, 



CO 



4* = p$+e>>, 

n=l 

A (n) = -J2^ kl VS k2 ■■■VS kn+1 , (29) 

(n) 
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where we introduced the notations (k > 1) 



qO _ p(0) qk 

o — r i N , o 



I -P. 



(o) 

AT / ) 



(30) 
(31) 



(0 



fcl ..Km — 



and the prime on the sum in (I3"T|) implies fci + . . . + k m — I. The square root operator in 
eq. (12 8 p is defined from its series expansion, 



p(0) p p(0) 



-1/2 



m=l 



(2m- 1)!! 
(2m)!! 



n=l 



T() 



(32) 



where A*™) = PfQA^ 'P^°2- Up to and including fourth order in V we find 



p(0) p p 



(0) 



-1/2 



if) _ 1^(2) _ i^(3) + 

j.tv 2 2 



1 



3 ^(2)) 2 _I^(4) 

8 2 



(33) 



and corrections are of fifth order in the perturbation V . 

Likewise, the Kato-Takahashi operator (|28|) can be obtained in a series expansion. Up to 
and including the third order in the perturbation V, the Kato-Takahashi operator = 
E^°=o r£v reads 



f(o) _ 6(0) 

1 i,N — r i,N > 



<1) 
i,N 

,(2) 



SVR 



(o) 

,N ' 



1 



r# = -s 2 vP^vpw - ~P^vs*vpw + svsvp^ , 



,(o) 



,(3) 



>(0) 



(o) T >6(o) 



>(0)t>£(°)i 



,JV ~~ ^ i,N i,N i,N ' 2 i,N i,N i,N 2 ''^ ''^ 



,(0) 



- &vp!®vsvp!% - svs 2 vP^ T vPP 1 ^ (0) 



(34) 



,(o) 



iJV' 



i,N 



/..V 1 ' /. Y 2 ' 



1 x 

2 



~svp§vs 2 vp§ - -P^%vsvs 2 vP^ + SVSVSVP^ 



i,N ■ 



With these definitions we can express the eigenstates of the Hamiltonian H, \^f) 6 £o,l ; at 
half band-filling (N = L) in terms of the eigenstates of the Hamiltonian Hq, |$) e £o°2, 



l*> = r 0l r|$) • 



(35) 



We introduce the state 



7!' 



r) + 



(36) 
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which is a symmetric mixture of the two ground states of the single-impurity model in the 
atomic limit at half band-filling, 

(L-3)/2 

\<h) = d+ JJ ajajjvac) , (37) 
1=0 

(L-3)/2 
1=0 

Note that the overall phase of the states is fixed by the convention that an electron on a 
particular site with spin "j" is placed to the left of an electron with spin J,. Then, eq. ( 121 j) 
reduces to 

Glhb(w < 0) = ($|f+ L rf+(u; + # - £o,l - i^)- 1 ^ 0>i |$) , (39) 

which is an exact expression. The average can now be taken in the known ground states 
of Ho- The energy Eq^ belongs to the symmetric single- impurity Anderson model and 
should not be confused with the ground-state energy E (L) of the Hubbard model at half 
band-filling. 

2. Transformation of the multi- chain Hamiltonian 

In order to make use of eq. ( 1391 . we must project H onto the eigenspaces of H Q . Before 
we can continue, we must be aware of the fact that the lower Hubbard band consists of the 
primary lower Hubbard sub-band, centered at (—{7/2), and higher-order sub-bands, centered 
at (—{7/2 — iU) (z = 1,2,.. .). We assume that, (i), these bands do not overlap in the Mott- 
Hubbard insulator, and that, (ii), the degeneracies of the eigenspaces Ef^ are not lifted to 
all orders in perturbation theory, as is the case for the ground state at half band-filling. 6 
Under these assumptions, the operators P^l-i project onto those states, forming the ith. 
lower sub-band. Since these projectors form a complete set, J2i^i,L-i — 1) we ma y simplify 
the Green function (139]) as follows: 

= ^($|f+ L 4A,L-i(^ + H- E 0tL - ir ? )- 1 P iii _ 1 4fo, L |$) , (40) 
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where we used the fact that the Hubbard sub-bands do not overlap. Now that ~tf N t iy N = P^ 



and r^jvl^jy- = Pi ; N hold, we may further simplify (HUj) 

Glhb(^) = ^(^If+^f i , L _ 1 f+_ 1 (w + H- E , L - «7) _1 f yi-ifjudrf o, L |$> 

= J>|d+ (^-l + ^ - Vjl - ^"^tl^ ( 41 ) 
i 

with the 'reduced' operators 

^ = rJ x _ 1 d <T fo,i , = (di,a) + , (42) 
= f t L ^Hf hL ^ . (43) 

As seen from eq. (l4"Ii) . we need to work with the 'reduced Hamilton operator' h^L-i which 
describes the dynamics of a hole in the symmetric Anderson model. 

Kato's perturbation theory^ also provides a perturbative expression for the projected 
Hamiltonian, 

oo 

(h-E^P^) = (44) 



n=l 



£(n) ._ S kl VS k2 V . . . S kn VPl% . (45) 

(n-l) 

In order to evaluate (j4"3"]l we have to multiply the expression (1141) with (P^P^P^) -1 / 2 
from the right and with T^y from the left. The various orders in the interaction are then 
combined to give 

oo 

rt, N {H-E$)t itN = J2xU ( 46 ) 

n=0 

where, up to and including the third order in 1/U, we find 

6(0) _ 6(0) t7 6(o) 
i,N ~ r i,N v r i,N 

5(i) _ 6(0)t/ct>6(0) 

n i,N — r i,N v ° v r i,N ) 



(47) 



3$ = ^y£y£y£y^2 + ^ 



p( )t>c3t>6(°)t>6(o)t>6(o) , p(o)t7 6(o)t7 6(o)t7c.3t7 6(o) 



-P^VS 2 VSVP^VP^ - P^vs 2 vP^vsvP^ - P^vsvs 2 vP^vP^ 
-P^vP^vs 2 vsvP^ - P^vsvP^vs 2 vP^ - P^vP^vsvs 2 vP^ 
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From the operators R\ ^ we readily obtain the correction to the ground-state energy at half 
band-filling (N = L), 

oo 
n=0 

<t = (48) 

Note that Eq~^ = —(L — 1)17/ 2 — U/A because the [L — l)/2 sites in the lower chain are 
doubly occupied in \(f>^), the impurity site is singly occupied, and the upper chain is empty, 
cf. eqs. (1241 and (1371). 

In order to single out the leading-order contribution in eq. ( l4Tj) we use e\ — E^~^ = 
U/2 + iU and define 

oo 

4l-i = kL-i - (Eo,l + U/2 + iU) P^_ x = (49) 

n=0 

with 

^l-, = ^li - <^lx ■ (50) 

Then, the Green function for the lower Hubbard band fj4T|) can be written as the sum over 
the contributions from the individual sub-bands around Ui = —U/2 — iU, 

GlhbM = ^($|<t (V + U/2 + iE0^2_! + Ci,L-i - vq) ^ lA®} ■ (51) 

i 

For the primary Hubbard sub-band, we will drop the index i = whenever possible. 
B. Matrix representation of the self-consistency equation 

The self-consistency equation ( 1201) must be solved for all frequencies in the respective 
sub-bands of the lower Hubbard band. The density of states of the individual sub-bands 
can be viewed as probability distributions which can be characterized by their moments. 
The idea is to express the Green function and the hybridization function by their moments 
so that the self-consistency equation reduces to the condition that the two sets of moments 
agree. In this way, only a countable set of numbers must be compared. A suitable way to 
generate moments from a Green function is provided by the Lanczos iteration procedure. 
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1. Lanczos iteration 



For the starting vector |0) and the Hermitian operator O, we use the following form of 
the Lanczos iteration: 

|1> = -O|0) + ao|0) , 
\n+ 1) = -0\n) + a n \n) +b n _ l \n- 1) , n > 1 , (52) 



(n\ 


o\ 


n) 


(n\ 


n) 



a n ■= v 1 1 , n > 0, (53) 



where 



(n - l\0\n) (n\n) ^ 1 , A . 

bn-i := t ' \ = -7 \ 1N , n > 1 . 54 

(n — — 1) (n — l|n— 1) 

The matrix representation D of O within the Lanczos basis {\n)} is tridiagonal, symmetric, 

and we have 



On,n = ~i 1 \ = a n 5 (55) 



(n\ 


o\ 


n) 


(n\ 


n) 



(n — l\0\n) i — — . nS 

On-l,n = \ ' = " \f-bn ■ (56) 

y — l\n — l)y (n|n) 

In the following we use fracture letters for the matrix representations of the corresponding 
operators in their Lanczos basis. Note that the parameters b n can only be defined up to an 
arbitrary phase, which is a sign factor for real matrix elements. Thus, the matrix D', where 
we change the sign of an arbitrary off-diagonal element and of its symmetric counterpart, 
represents the same operator O. 



2. Hybridization function 

We introduce electron baths for every sub-band of the lower Hubbard band. Starting from 
the single- impurity Anderson model in star geometry, we write the hybridization function 
in the form 

oo t^2 oo 

where Aj(w) denotes the contribution of the ith sub-band. We cast each A, (a;) into matrix 
form by applying the Lanczos iteration with the starting vector 

y gi m 
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where gi is the weight of the ith sub-band in the density of states, J2i9i = 1/2, see ap- 
pendix |X] With the (discretized) Hamiltonian for the bath electrons in star geometry 



(59) 



we may write 



-i 



Ai(u) = i(0\ {(u + U/2 + iU) - H A>i - i V ) |0), 
= {{{u + U/2 + iU)t-^ i -ir,y 1 ) m . 



(60) 
(61) 



This form can be verified by noting that [HA,i] n \0)i = 1/ \Z^^^^tji 1 c7 ^ // ^)^i(^*2;^^)^^'^77i ; crl^^'^)* 

We note that the mapping of the single-impurity model from the star geometry to the 
multi-chain geometry is based on the Lanczos procedure. 3 Therefore, the starting vector |0)j 
is identical to an electron at the first site, |0)j = (l/\/2)a£ J vac), of the ith lower chain 
in the multi-chain geometry. Thus, the Hamiltonian H A i for the bath electrons can also be 
written in the form 



1=0 cr 

Then, the matrix t)A,i representing Ha% in the Lanczos basis reads 



(62) 



/ 



\ 



(63) 



£i;0 U;0 
ti;0 ti-l 

ti;l £i;2 ti-2 

V '• / 

where the entries not shown are zero. The parameters e^ m and t^ m define the single-impurity 
Anderson model in its multi-chain geometry. 



3. Green function 



For the Green function (ED we use 



as the starting vector, see eq. (|36j) . and 



|*o) = d t |$) 

6 = A,l-i 
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(64) 



(65) 



as the operator in the Lanczos iteration, see eq. (150|) . In this way we obtain the matrix 
representation 

oo 



' 00 



i=0 



where the structure of £; is given by 

/ 



Ci;0 T i;0 

n-o e i; i n-i 

Ti-1 Ci-2 Ti-o 



\ 



\ 



J 



The parameters e i]m and r i;m must be calculated from eqs. (157)1) and (1B^|) . 



(66) 



(67) 



^. Self- consistency equation 

For t = 1 as our unit of energy, the self-consistency equation fliZOl reads (u; < 0) 

oo oo 

(((" + + *E7) 1 - t) Aii - ir?)- 1 )^ = Y ((( w + f 7 / 2 + «0 * + £ ~ ir/)^) 00 - (68) 

In this work we are mainly interested in the primary lower Hubbard band. As we show in 
appendix [Aj it is the only lower Hubbard band with non- vanishing spectral weight up to 
and including third order in 1/U. Up to this order, we may therefore write 

(((w + U/2) t-U- ^) _1 ) 00 = (((w + U/2) 1 + £ - i^)- 1 )^ , (69) 

where we have dropped the subscript i = 0. Reckoning with (1691 . we realize that 

f)A = -£ (70) 

is a sufficient condition to ensure the self-consistency. From the continued fraction expansion 
of the hybridization function and the Green function it can readily be shown that it also is 
a necessary condition. Therefore, the self-consistency condition reduces to 



-n &n and |t n | 



(71) 



We already remarked that the off-diagonal Lanczos parameters are only defined up to a 
sign factor. Note that eq. ( ITTj) is a vast simplification over (T2UI) because, due to the matrix 
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structure, we only have to equate numbers and not functions. However, since all our calcu- 
lations are implicitly done in the thermodynamic limit, there is a (countably) infinite set of 
parameters to fix. As we shall show explicitly up to third order in 1/U, the locality of the 
Hubbard interaction guarantees that there is an index l m in mth order perturbation theory 
such that the Lanczos parameters r n and e n become constant for n > l m . 



IV. SOLUTION OF THE DMFT EQUATION 

In the first part of this section we show how the DMFT equation is solved to leading 
order in 1/U. In the second part we summarize the results to third order. 



A. Calculations to leading order 

First, we calculate the ground-state energy and set up the starting vector for the Lanczos 
iteration. Next, we calculate the action of the Lanczos operator on the states with a single 
hole in the lower Hubbard chain. Then, we derive the parameters e^ °\ s^, t °\ and if 
from the first two Lanczos iterations. Lastly, we prove the result = and t\ = 1 for 
all I by induction. 



1. Ground- state energy and starting vector for the Lanczos iteration 

To leading order, we have = from ( I34p . From f[3"o"j) we see that the ground state 
at half band-filling is not transformed to leading order, |\I/ < - ' ) ) = |<3>), see eqs. ( 136]) - (138]) . The 
correction to the ground-state energy to leading order follows from eq. (jlSjl as 

(L-3)/2 

4 o L = (0ti4 o) i0t) = (0ti^i0t) = (0ti^i0t) = 2 £ i- ( 72 ) 

1=0 

The operators V and V\ do not contribute because they modify \4>\)- The contribution of V 2 
is readily calculated because the sites of the lower chain are doubly occupied, see eq. (1251) . 
According to (|64p . the starting vector for the Lanczos iteration is given by 

rt 0) ) = p1 o) A^ o) I$)-I0-i), 
/T- at( l - 3 )/ 2 

i0-i> = y 2^ti0t) = y o n & tt & tj vac ) ■ ( ?3 ) 

V V 1=0 
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Note that, in general, the starting vector is not normalized to unity. 



2. Lanczos operator 

The operator for the Lanczos iteration (165]) is given by 

41 = - E«m, = (l> - <L) Pfl, , (74) 

see eqs. (jUJ) and (loll . In the course of the calculations, we shall need the eigenbasis of P^ln 
i.e., single-hole states in the half-filled ground states of Hq. Apart from in ( 1731) . we 

define for n > (see eqs. f )37|) and (138]) ) 

(L-3)/2 



:= -J~dfa n A ] ] d ^ | vac) , (75) 



2 

(L-3)/2 



Z=0 

and their useful linear combinations 

l7n) := (-l) n y|(|0n; U )-|0n;,)) , 

:= (-l) n y|(l0n;«) + IXn)) , (76) 
K;rf) - (-l) n yf (l0n; d > + |Xn)) • 

Note that the states are not normalized but they are site-orthogonal in the sense that the 
overlap between states with different site indices is zero. 

The action of the Lanczos operator on the states is readily calculated. One finds 

for n = 0, 

41i0-i> = Ito) , 

4°lil7o) = |0-i)-4 O) |7o) + 4 O) |7i), 

4°-ll m 0;«) = 2^-l) _ £ 0) | W 0;«) +4 I m l;«) > 

4°-il m o ; d) = — — [0_i> -e[, 0) |mo;d) +4 0) |mi ;d ) , (77) 
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and, for n > 1 and x n = j n , m n]U , m n;d , 



The effect of the Lanczos operator is identical for all n > 1. Note that this holds true in 
mth-order perturbation theory in 1/U for n > m + 1. 

3. First and second Lanczos iterations 

In the first Lanczos iteration we must determine the state 

1^):=-^^) + ^^) (79) 

with 

" <*n*r> ' ( ' 

With the help of (E7J) we find 

lM 0) ) = -|7o) , ej 0) = (81) 

because = 2(4>-i\lo) = 0. The self-consistency equation (j7T|) then gives eg ' 1 = 0. 
Furthermore, we have 

° = (ifiif) = <70 ' 70> = ' ' ' 

so that = 1 follows from the self-consistency equation (17X1) . 

In the second iteration we can use the results from the first iteration. We calculate 

|*f ) := -All*?)} + eS 0) |M 0) ) + r ( ?W> (83) 



2«7o|0-i) + <7o|7i» = O, (84) 



with 

(0)|/S(0) ,^(0 v 
(0) l^-L-ll^l / 

61 ~ <*i 0) l*! 0) > 

where we used £q = and = 1 in (ITT)) . The self-consistency equation (J7TJ then gives 
= 0. 

The second state in the Lanczos iteration reduces to 

l^ 0) ) = l7i> , (85) 
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and we find 

(0) \ W 1 ^L-l W 2 / 



so that = 1 follows from the self-consistency equation 



2((<M7i> + <7i|7i» = -l, (86) 



4- Induction 

Now we are in the position to prove by induction that e„ = = and Tn^ = — 1 = 
— tffl- Let this assumption be true for 1 < n < M — 2 (M > 3) and assume for 1 < n < M — 1 
(M > 3) that 

l*! 0) > = ("l) n |7n-i) • (87) 

Then, we calculate 

I vl>< )\ — lvl> (0) \_Up (0) IvI/ (0) \_Ut (0) IvI/ (0) \ 

with 

/^(O) 1/5(0) ,^(0) v 

eV-i = — ( 0) (0) — = 2(-l) + ((7M-2I7M-3) + (7m-2|7m-i)) = , (89) 

where we used s^_ 2 = and t$_ 2 = tf^_ 3 = 1 in (!7H|) . The self-consistency equation ([71]) 
then gives e^_i = which proves the next step in the induction for 6n\ 
The next state in the Lanczos iteration reduces to 

I*a?) = (-1) M (i7M- 3 ) + |7a/-i» - (-l) M - 2 |7A/- 3 ) = (-1) M |7m-i> , (90) 

which proves the next step in the induction for |\l/i ^). 
Finally, we find 

/viv(0) 1/5(0) I o,(0) v 

r M-i = (o) (o) = ~2(-l) «7Af-3|7Af-i> + (7A/-i|7Af-i)) = -1 , (91) 

so that t^-i = 1 follows from the self-consistency equation ( ITT]) . This proves the next step 
in the induction for t^P ■ 

We recall that our approach strongly relies on the simple form ( 120]) of the self-consistency 
equation. For a general form of the bare density of states p(u>), the calculation of the density 
of states for the lower Hubbard band to leading order is a demanding task.— 
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B. Results up to third order 



The calculations up to third order are straightforward but tedious. Simplifications arise 
from the fact that we are interested in the half-filled case (N = L) and the situation with 
a single hole (N = L — 1). Moreover, the results to leading order simplify the analysis 



considerably. Details are given in Ref. 20] . Here we summarize the results. 



1. Ground- state energy and starting vector for the Lanczos iteration 

The ground-state energy is given by 

U 1 3 (L_3)/2 / 1 \ 

Eoj. + U(L-l) + - = (<h\R L \<h) = ----=+ £l + [m) ■ ^ 

The starting vector for the Lanczos iteration to third order is given by 

l^o) = - -jj\m . u ) + + \m 1;u )j - ^ \^\m . u ) + |m 2;u )^ . (93) 

In the actual derivation, the results of the lowest-order calculations are used in first order, 
those of the first-order calculations are required in second order, and so on. 

2. Lanczos operator 

Up to and including the third order in 1/U the Lanczos operator Cl-i = £ reads 

(L-3)/2 



V 1=0 ' 



~h + ^ (k ~ \kk - h h^j (94) 
+ ^3 (~^ 3 + h2h ° + hofl2 ~ \ h ^ho? ~ ^(ho)% + (hi) 2 



Here, we used the abbreviations 

>(0) T 7 6(0) 



h — p « t/ pw 

h = P^VoSVoP^ , 
h = P^VoSVoSVoP^ , 

h = P^VoSVoSVoSVoP^ , (95) 
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where Vq — Vq + V^ describes the coupling of the chains to the impurity with amplitude 
Vq = l/\/2 and the free hole motion along the chain with amplitude t\ = — 1. The operator 
S measures the inverse number of excitations above the ground states of H for N = L — 1 
particles, 

oo p(0) 

i=i 3 

The remaining task is the calculation of the action of the Lanczos operator on the single-hole 
states fl75|) and fl76|) . 

Up to and including third order in 1/U we find for the state with the hole at the impurity 

£|« = { 1 -ip) w + (h + w) ^ + ^ w ■ (97) 

For the states with the hole at site n we find for n = 

1 |T2> , (98) 



AU 3 

2 JO 2 (J) 



1 / 3 \ l l £ w £ w 

^ ' 1=0 1=0 

1 / 3 \ i i 2 e O 2 



=0 

For n > 1 we find 

11 1 



+ E ^7 (*n-ll7n-!> - eW| 7 n> + fJWO) , 

1 2 1 

>C|m n;u ) = -5 nA -^\m Q . d ) (*i-l|»"n-l,ii) - 4°l m n.;«) + 4 Q | m n+l;«)J , (99) 



«=0 
2 



t\m n - d ) = S n)1 -—-^\m 0lu ) + E 771 (*n-ll m »- 1 ;<*) _ 4^l m n;d) + t { n\ m n+l-,d)) ■ 

1=0 

In the above formulae, we did not include the third-order contributions to \m n - u ) and \m n -d) 
because they are not needed for the third-order calculations. 

With the help of the starting vector (l93"j) and the action of the Lanczos operator f l97l) - 
( 199|) . the Lanczos vectors can be generated iteratively, along with the values for e; and r/, 
see eqs. ([55]), (121, and flBT]). 
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3. First and second Lanczos iterations 



The first and second Lanczos iterations give the vectors 

= — |To> + jj [~\^-\) + K;«>) + Jp (jl7o> - ^\m ; U ) ~ \m 2 -u) 
1/9 3 25 \ 

+ Jp " I l7l) + Y |mi;u) + |m3;u) J ' (100) 

l* 2 ) = |7i> + \j Ql7o) " K«>) + ^ (-^-1) - |lTi> + \\mv, u ) + \m 3 Sj 

1 /ll, . 1. . 1. . 7. A . . 

+777 7rl7o) + oN) - >0;«) - d™2 ; «) - |m 4; „) • ( 101 ) 



f/ 3 V 4 " 7 2" 1 T u ' 2 
Moreover, the diagonal and off-diagonal Lanczos parameters become 

(*o|£|*o> 7 



e = 



(#o|*o> " 4C/ 3 ' 



ei = 



1 31 



217 8C/ 3 ' 



ON 


*2> 


<*ll 





For n > 3 the Lanczos matrix elements become constant, and the Lanczos vectors obey a 
building principle. This can be proven by induction. 



4- Induction formulae 



Up to and including third order in \jU we have 



en = e» = -^-J, forn>2, 

-t„ = r n = -l-^ forn>l, (104) 



and, for n > 3, 



(-l) n |^ n ) = | 7 n-l) + ^ Q|7n-2> - \m n]l 



+ (^—lln-s) + OnlTn-l) + ^|™n-l;u) + |™ n+ i ;u )^ 

1/1 11 

+ ^ ( g |7n-4) + &n|7n-2) + ^\ln) ~ ^\m n ^ 2 ;u) - C n \m n - U ) - \m n+2]l 



(105) 
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with a n = 3(n 
|T-i> = 



3)/8, b n = 3(n - 4)/16 + 27/8, and c n = 3(n - 4)/8 + 17/4, where we set 



The induction proof is lengthy but straightforward and can be found in Ref. [20 ]. 



V. HUBBARD BANDS IN THIRD ORDER 



In the last section we calculated the parameters for the hybridization function of the 
single- impurity Anderson model in two-chain geometry up to and including third order 
in 1/U. The DMFT self-consistency equation on the Bethe lattice (12"U|) shows that the 
hybridization function is identical to the impurity Green function which, in turn, is equivalent 
to the local Green function of the Hubbard model. 



A. Density of states of the lower Hubbard band 

The hybridization function A(u) can be obtained from an equivalent scattering problem 
for a single particle on a semi-infinite chain. We start this section by formulating this 
problem. Next, we calculate the single-particle gap and the hybridization function in closed 
form. Lastly, we expand this expression systematically in 1/U which defines the 'band-part' 
Green function. 



1. Scattering problem 

The Lanczos algorithm provides the tridiagonal matrix representation of the hybridization 
function, 

/ 

t S\ t 

let 



\ 



V 



(106) 



where 



to = 1 + 



AU 3 
1 

81P 



1 31 

+ 



1 



35 



' £l 2U ' £n 2U 8U 3 



= e (n > 2) 



(107) 
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Note that only odd (even) orders appear in the l/?7-expansion of E\ {t\). 
As shown in Sect. IIII B\ the hybridization function can be obtained from 

AlhbM = (((w + U/2) l-t)A~ ^) _1 ) 00 • (108) 

The matrix 1}a H106[) corresponds to a tight-binding Hamiltonian K which describes the 
transfer of a single particle on a semi- infinite chain, compare eq. (I25I) . plus a scattering 
potential W at the boundary of the chain, 

H scat = k + w , 

oo oo 

k = tj2&)(i+i\ + \i+i)(i\)+£j2\ l )( l \> ( 109 ) 

1=0 1=0 

w = ^|o)(o| + ^|i)(i|+tS(|o)(i| + |i)(o|) 

with e* = e - s = -1/(2U) - 21/(8f/ 3 ), ej = e 1 - e = -1/(2U 3 ), t* = t - i = -1/(AU 2 ). 
Note that the scattering potential W is attractive which results in a redshift of the density 
of states; see below. 

The hybridization function can equally be calculated from the one-particle Hamiltonian 

AlhbH = (((u; + [//2)l-f) A -i7 / )- 1 ) 00 

= 1(0| ((w + f//2)l-4ca t -ir/)^|0) . (110) 

In this way, the Green function for the lower Hubbard band G ? lhb(w) can be deduced from 
a one-particle problem. 



2. Single-particle gap 

The attractive W is too weak to generate a bound state below the lower band edge of 
the tight-binding operator K. Therefore, it does not change the support of the imaginary 
part of the bare Green function defined by K which is given by \u — e\ < 2t. In turn, this 
implies that the upper edge for the lower Hubbard band is given by = — U/2 + e + 2i so 
that the charge gap in (|22|) is given by (bandwidth W = 4) 

A c = 2|- C //2 + f - + 2f|= C ;-4-I-^-| J + (i I ) . (Ill) 
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The result to second order was derived earlier by Eastwood et alii.S Note that the coefficient 

n 

to third order is larger than anticipated in Ref. [6]. 

Let A C (U C ) denote the critical value of the on-site interaction U where the charge gap 
closes, i.e., A C (U C ) = 0. Up to third order we find 

UP = 4 , 

£/ c (1) = 4.236 [5.90%] , 

U^ 2) = 4.313 [1.82%] , 

U^ 3) = 4.406 [2.16%] , (112) 

where the number in square brackets gives the percentage change to the result of the previous 
order. Apparently, the changes in the estimated critical interaction strength from the second 
to the third order are of the same order of magnitude and the critical interaction strength 
does not converge quickly in these low orders. 

3. Green function for the scattering problem 

The calculation of the boundary Green function for a semi-infinite chain with a local 
potential at the boundary is readily accomplished.— 1 ^ We define the general two-site Green 
functions 

g l>m (cu) = (l\ (( W + U/2) 1 - # scat - iTj) _1 \m) , (113) 

s£M = CI (V + I 7 / 2 ) 1-K-ivy 1 \m) • (114) 

The Green functions (II 14ft for the tight-binding Hamiltonian K are calculated explicitly in 
appendix [B] 

In ( TTT3]) we use the operator identify (A - = A' 1 + A~ 1 B(A - B)' 1 with A = 

(co + U/2)l — K and B = W so that we can write 

oo 

So,oM = 4oM + E 9$(u)(l\W\m)g m>0 (u) (115) 

l,m=0 

= 9o°fi{u) + erfoM^o.oM + elgj$(u})gi fi (u)) + t* (g§ ] Q (uj)g lfi (uj) + #JiM#o,oMj , 
where we used the locality of the scattering potential W (11091) in the second step. Likewise 
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we obtain 



9i,o( u ) = #i°d( w ) + e o0i"o( w )0b,o(w) + e^i"i(w)pi,o(w) + t* ( ^(u^o.oH + 0i"o(w)0i,o(w 



,(0) 



(116) 



We can solve the coupled equations HI 1 5j) and (11161) to give our final result (g 00 (uj) 
2G LHB M) 

0o?o( w )0i°iM - ^o°i ) ( w )^°d( w ) 



s8(w)-£i 



#0,0 J 



iV(w) 



+ (eSe! - (^) 2 ) (<dM<??iV) 
From appendix [B] it follows that gf'l (oj) = g^l (oj) and 



(117) 



(118) 



so that we can cast our final third-order result into the form 

n 2 



0o,o(w) 



5 + ej) g$(u) + (<■•<■• - (t*) 2 - 2t*t) 9 ™( u ) 



*+2 



e\t 



3 ■ 



(119) 



The density of states is the imaginary part of this expression, 27iD L h B (uj) = lm[g 00 (u)}. The 
bare boundary Green function is given by [x — (u + U/2 — e)J (2t)] 



tg$(u}) = B(x 2 - 1) (x - sgn(x)Vx 2 - l) + (l - x 2 ) x + iVl 



(120) 



where Q(x) is the Heaviside step-function. For the density of states we only need the region 
\x\ < 1. 

In Fig. [2] we show the results for the density of states of the lower Hubbard band for 
U — 5 (bandwidth W — 4) to first, second, and third order in 1/U. The overall spectra 
display a redshift of the Hubbard semi-ellipse f fTOj) which describes the density of states to 
leading-order, 2D^ b (uj) = p(co + U/2). The spectra to higher orders differ from each other 
mostly by a shift in the spectral support so that the deviations are best visible close to the 
band edges. 
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FIG. 2: Density of states of the lower Hubbard band, ttD^ b (oj) for {7 = 5 (bandwidth W = 4) 
up to and including orders n = 1, 2, 3 (black, blue, red colors). 

4- Band part of the Green function 



The full solution ( I117P contains higher-order corrections in 1/U due to the interaction- 
dependence of the denominator N(u). We may expand it order- by-order to derive a Taylor 
series in 1/U for the Green function. Such an order-by-order expansion ignores the fact 
that the attractive potential W generates resonance-contributions at the band edges of the 
Hubbard band; see below. Therefore, we denote the Green function from the order-by-order 
expansion as 'band-part' Green function. It can be cast into the form 

3 



n=0 



2 1 + 



X 



00 + 



U 



1 



35 



8U 2 J 2 2U 8U 3 ' 

g n (x) = 6 (x 2 - l) (r n+ i(x) - sgn(x)Vx 2 - lU n (xj) , 

9n (x) = Q(l-x 2 ) (T n+1 (x) + iVT^U n (x)^) , 

with T n (x) [U n (x)} as the Chebyshev polynomials of the first [second] kind^ 3 . and 

1 39 



121) 



A, 



— 1 > Ai 
1 



4U 2 



A, 



2U 1QU 3 
1 

'8LP 



(122) 
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are the expansion coefficients. The first-order result was derived earlier in Ref. [28|. Using 
an intuitive method, Eastwood et al.- derived the 'band-part' Green function to second order 
in 1/U for the Hubbard model in infinite dimensions. So far, their method could not be 
extended systematically to higher orders. 
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(a) n — 1 
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(c) n = 3 



FIG. 3: Resonance contribution to the density of states of the lower Hubbard band, D^ b (uj) = 
-Dlhb(^) — D^^(uj), as a function of frequency in nth order perturbation theory for U = 5.5 
(bandwidth W = 4). 



In Fig. [3] we show the resonance contribution to the density of states, D™^ b (uj) for U 



5.5. 



It is defined as the difference between the band part ( 112111 . D^ B (u), and the full density 
of states -Dlhb(cl') ( II 19fl . The difference is seen to be fairly small which had to be expected 
because the potential W is rather weak. In general, the resonance contributions increase 
slightly the density of states close to the band edges and decrease it in the middle of the 
band. 



B. Comparison with numerical results 

Finally, we compare our analytical results with data of advanced numerical methods 
for the DMFT for the Mott-Hubbard insulator. The Dynamical Density-Matrix Renor- 
malization Group (DDMRG) method provides the gap and the density of states at zero 
temperature. 8 Quantum Monte-Carlo (QMC) gives the Matsubara Green function at low 
but finite temperatures. Ref. (| contains a comparison with early methods in the field. 
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FIG. 4: Charge gap as a function of the interaction strength for various orders in the l/C/-expansion. 
The dots are DDMRG data points.^ 

1. Gap 

In Fig. H] we show the g function of the interaction strength for various orders 



in the 1 /[/-expansion together with the DDMRG data of Ref. [8J. The third-order theory 
reproduces the DDMRG data points very well. 

Note, however, that in another DDMRG study^ the gap closes around U = 4.8. The 
differences in the two approaches lies in the reconstruction of the density of states and 
the extrapolation of the gap from the finite-size data. Apparently, different reconstruction 
algorithms can result in substantially different extrapolations close to the transition. 



2. Lower Hubbard band 



In Fig. Owe show the density of states for U = 4.8 to third order in 1/U together with 
the DDMRG data of Ref. [gj]. The overall agreement is very good. This has already been 
observed from the results to second order.- 

It is seen, though, that a resonance develops at the upper band edge in the DDMRG 
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(a) U = 5 



(b) U = 4.8 



FIG. 5: Density of states of the lower Hubbard band from third-order perturbation theory in 1/U 
(full line) in comparison with DDMRG data points 8 for U = 5 and U = 4.8. 

data which is not seen in perturbation theory to third order. For U = 4.5, the resonance 
is more pronounced^ and resembles the split quasi-particle peak of the metallic phase. One 
may wonder whether such a resonance could be obtained from higher-order perturbation 
theory. A model study 32 shows that the parameter set = —0.2 and £l< m <g = 0.1/m in 
the scattering potential W can readily account for both the overall redshift of the density of 
states and a resonance at the upper band edge. Since the range of the repulsive potential is 
finite, an expansion of the density of states to high but finite order could possibly reproduce 
the resonance seen in the DDMRG data. 

3. Matsubara Green function 

The Matsubara Green function for the Hubbard model is defined by 



where /3 = l/k^T is the inverse temperature, Q is the grand-canonical potential, and T T 
orders the operators in imaginary time. The operators in imaginary-time Heisenberg repre- 
sentation are defined by (—(3 < r < j3) 




1 



(123) 





(124) 



32 



-14 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 



T 



FIG. 6: Third-order result and QMC data for the Matsubara Green function for U = 6 (blue) and 
U = 5.2 (green). The inverse temperature is j3 = 20 (T = 0.05), the gaps are A C (U = 6) = 1.75 
and A C (U = 5.2) = 0.89, respectively. Note that the data are shown on a logarithmic scale. The 
shading indicates the statistical error in the QMC data. 

The Fourier transformation of the Matsubara Green function is denned on the points iu n = 
(2n + l)n/j3 (n: integer) on the imaginary axis. 

The retarded Green function at finite temperature T is obtained from the analytic con- 
tinuation, 



Therefore, we may express the Matsubara Green function with the help of the density of 
states at finite temperature in the form 



with < r < p. Note that it is easy to evaluate (I126P for a given density of states but it is 
very difficult to reconstruct the density of states from numerical data for Q{r). 

For very low temperatures and for large interaction strengths we approximate the density 
of states by its zero-temperature expression to third order, 



a 



' Iet (u; (3) = Q{iuj n ->• u + ir]) . 



(125) 




(126) 



oo 



e 



du .Dlhb(w) + D LllB (-u) 



. e-P" + 1 



(127) 



— oo 
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FIG. 7: Difference between the third-order result and QMC data for the Matsubara Green function 
for U = 6 (blue) and U = 5.2 (green). The inverse temperature is /3 = 20 (T = 0.05), the gaps are 
A C (U = 6) = 1.75 and A C (C7 = 5.2) = 0.89, respectively. Note that the difference is augmented by 
a factor of 10 3 to make it visible. The shading indicates the statistical error in the QMC data. 

which we compare with QMC data of N. Blumer.— The approximation is not as drastic as 
it may seem because, deep in the Mott-Hubbard insulator, thermal excitations are expo- 
nentially suppressed due to the finite charge gap. Therefore, corrections to f!127j) should be 
exponentially small in A c (U)/k^T. 

In Figs. [6] and [7] we compare our analytical results (I127P to QMC data for (3 = 20 (T = 
0.05) at U = 6 and U = 5.2 where the gaps are A C {U = 6) = 1.75 and A C {U = 5.2) = 0.89, 
respectively. The results agree very well. Note, however, that Q{r) is rather feature-less so 
that fine points such as the width of the Hubbard bands or the density of states cannot be 
reconstructed easily from QMC data for Q(t). 

VI. CONCLUSIONS 

In this work, we have studied the Mott-Hubbard insulating phase of the Hubbard model 
on a Bethe lattice with infinite coordination number. We have adopted the Kato-Takahashi 
perturbation theory to solve the self-consistency equation of the Dynamical Mean-Field 
Theory for the symmetric single-impurity Anderson model in perturbation theory up to and 
including third order in the inverse coupling strength U . To this end it has been necessary 
to use the mapping of the single-impurity Anderson model from the 'star geometry' onto the 
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'two-chain geometry' which represents the energetically separated lower and upper Hubbard 
bands. In higher orders, a multi-chain mapping is required in order to resolve the various 
Hubbard sub-bands. For the present study, we could ignore the secondary Hubbard bands 
whose weight is of fourth order in 1/U . 

Our results for the Mott-Hubbard gap reproduce those of an earlier analytic study 6 - of 
the Hubbard model on the Bethe lattice with infinite coordination number. We confirm 
the second-order results of Ref. js] and extend them to third order systematically. The 
agreement between the perturbation theory in 1/U and the DDMRG data of Ref. js| for 
the gap is very good. Note, however, that the precise value of the critical interaction U c 
where the gap closes, and the analytical behavior of the gap as a function of U close to the 
transition are still under debate.- 

The previous study 6 provides the Green function as a Taylor expansion in 1/U whereas 
the present study includes resonance corrections. The full density of states results from 
the calculation of the boundary Green function for a particle on a semi-infinite chain with 
nearest-neighbor electron transfers and an attractive interaction at and near the boundary, 
whose parameters we derived to third order in 1/U. For all interaction strengths where 
perturbation theory is applicable, U > 5 (bandwidth: W = 4, 4.4 < U c < 4.8), the 
resonance contributions are small. 

For U > 5, the agreement between the analytical results for the density of states and 
the DDMRG data& is very good for all frequencies. In addition, our zero-temperature 
expressions for the density of states provides a very good approximation for the density 
of states at small but finite temperatures. This can be seen from the excellent agreement 
between our approximate Matsubara Green function and Quantum Monte-Carlo data.— 

As in all kinds of perturbation theories, the number of terms to be calculated rapidly 
increases with the index of the order. In principle, the fourth-order terms could still be 
calculated 'by hand'. This requires a four-chain geometry so that the secondary Hubbard 
sub-bands can be treated properly. To fourth order there are more than 30 terms in the 
Kato-Takahashi operator and in the projected Hamiltonian. According to our analysis, 
much higher orders are needed to reproduce a resonance feature seen in the DDMRG data£ 
at the upper band edge of the lower Hubbard band. Such high-order calculations for the 
density of states appear to be forbiddingly costly within the DMFT. 

The ground-state energy of the Hubbard model on the Bethe lattice with infinite coor- 
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dination number was calculated to high orders using a computer algorithm based on the 
Kato-Takahashi expansion.— In the future, we plan to devise a similar algorithm for the 
calculation of the Mott-Hubbard gap. With a high-order expansion for the Mott-Hubbard 
gap we should be able to locate U c with a much better accuracy. 
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Appendix A: Weight of the secondary lower Hubbard band 

We apply particle-hole symmetry and the self-consistency equation (I20I) to the sum rule 
for the density of states^ and find 

i i t-0 00 r -1 00 

o = - / dudm [A(w)] = £ YsKn = E ^ . ( A1 ) 

^ " J — CO .■ n - _ ., -• n 



■- m 



i=0 



i=0 

where we use eq. ( IBTl) in the second step. The ith sub-band of the Hubbard band contributes 
the weight 

The weights gi can be calculated perturbatively. From the definition of the Green function 
of the lower Hubbard band (139|) we can readily write the sum rule (lAlj) as 

\ = - f (Mm [G L hb(w)] = ($ |f + L d+d t f ,l|$> • (A2) 

The state |\&) = d^Y q,l\&) can readily be calculated from the series expansion of the operator 
r ,L, eq. ( ESJ) , applied to the state |<3>), see eqs. (|36l)-(l38l). Up to and including third order 
in 1/U, we find 

3SLil*> = (l - ^2) I'M - ^ (l + ^2) Ku> + - i^|m 2;u ) (A3) 

in the subspace of the ground-states of Hq for L — l particles on L lattice sites. In addition, 
there is a finite second-order contribution in the subspace with excitation energy U above 
the ground states of H , 

P^Lm = j~ 2 (\Wo-,u) + \\^ d ) + lxS>) + O (^) , (A4) 



36 



where 

r- (L-3)/2 

V^oVo,, n + 



l0O;n) = \/5#uAu 11 & l,t & lJ VaC ) 

1=0 

(L-3)/2 



\(f>* . d ) = -W-^+ t ao jt JJ a+ d+ |vac) , (A5) 

/I. ( L - 3 )/ 2 

lxo> = -y 2^ a °'t II "tt d ul vac ) • 

Therefore, the weight of the secondary sub-band, i — 1, is given by 

* = <f|ASLl*> = 2p(j + i + i) = Ji. < A6 > 

and corrections are of higher order in 1/Z7. Because the higher sub-bands are even smaller 
in weight, (?j> 2 = 0(U~ 5 ), we can conclude from the sum-rule that g = 1/2 — 3/(4U 4 ). 
The direct calculation of g from eq. (1A3[) is not possible because it lacks the fourth-order 
correction — 7/ (4t/ 4 ) |</>_i) . 

Equation (1A6I) shows that, up to and including third order in 1/U, only the primary 
Hubbard sub-band contributes to the density of states for uj < 0. 

Appendix B: Green functions for a semi-infinite tight-binding chain 

Let B = YlHo 10 + M + K + -0(^1 describe the motion of a single particle over a semi- 
infinite chain. By definition, we have B\0) = |1) and B\l) = \l + 1) + \l — 1) for Z > 1. By 
induction we can prove the following lemma: 

(i) For all n £ No we find 

U n (B/2)\0) = \n) . (Bl) 
Here, U(x) are the Chebyshev polynomials of the second kind.— 

(ii) For all n, I £ N we have 

n 

U n (B/2)\l) = Q2(i-n)+2i\l -n + 2i). (B2) 

t=0 

Here, 9/ denotes the discrete unit-step function 6 : Z — > {0, 1}, 

, 1 for / > , 
<df-={ - (B3) 

for I < . 
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For the imaginary part of the bare Green function (11141) we can write (u = u + U /2 — e) 



Im 



gfl{u) =n(l\6[ut-iB)\m) 



which is finite in the interval \u\ < 2t. We set x = u/(2t) and obtain 



Im 



7T 



&<») =^(l\8[xl-JB/2)\m) 



(B4) 



(B5) 



When we formally expand the 'function' f(x) = tt5(x — z) (\x\ < 1, \z\ < 1) in a Chebyshev 

2D 

series,— 

oo 

n S{x -z) = 2Vl-x 2 J2 U n{x)U n {z) , 

n=0 

we can write the imaginary part of the bare Green function as 

oo 

: ^9 (1 - a; 2 ) VT^Ys U ^ x )( l \ U n(B/2)\m) . 

n=0 

Because B is Hermitian so that gf^u) = g^^u), we can restrict ourselves to I — m + h 
with h G No. With the help of the lemma, it is not difficult to show that for h,m G No, 



Im 



(B6) 



(B7) 



tlm 



,(°) 

m+h,m 



UJ 



Q(l - x 2 )VT^{5 mfi U h (x) 

oo 

+ (1 - 8mfi) ^2 ®rn+y-k<dk-y [U 2 k(x)6h,2y + U 2 k+1 (x)5h,2y+l] } (B8) 



k,y=0 



with the abbreviation x — (to + U/2 — e)/(2i), the discrete unit-step function 6^, see flB3j) . 
and <d( the Heaviside step-function. 

The real part follows from the Kramers-Kronig transformation.— We find 2 ^ 



tRe 



,(o) 

m+h,m 



UJ 



8mfllh(%) + (1 — ^m,o) ®m+y-k®k-y [hk(x)5h,2y + hk+l( x )8h,2y+l] ■ 

k,y=0 

(B9) 



As is proven by induction in Ref. 



20|, we have (n G N ,p G N), 



In( x ) = T n+ i(x) for |x| < 1 



I n (x) = T n+ i(x) — sgn(x)Vx 2 — lU n (x) for \x\ > 1 , (BIO) 
[I n {x)} p = T p(n+ i)(x) -sgn^Vx 2 - lC/ p(n+ i)_i(x) for \x\ > 1 , 

where T n (x) are the Chebyshev polynomials of the first kind.— In particular, the bare bound- 
ary Green function reads (a; < 0) 

tg$(u) = 6 (x 2 - l) (x - sgn^Vx 2 "^!) + 6 (l - x 2 ) \x + iy/l - xA . (Bll) 
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Finally, we note that powers of the bare Green function obey for n e No,p G N 



Fg%{u) = Q(x 2 - l)P n {x) + 9(1 - x 2 ) \T p[n+1) {x) + iVT^U p{n+iyi (x) 



(B12) 



With these relations and the properties of the Chebyshev polynomials, one can readily prove 
eq. ( DUD ; for details, see Ref. [20 1. 
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